      ##MMAD-RA##

#13 10 22

#Syntax of tables and figures in 'Introducing the MMAD Repressive Actors Dataset'

#Content:

#Load packages
#Import data
#Table 1 (page 5) - omitted
#Figure 3 (page 6)
#Figure 4 (page 7)
#Figure 5a & b (page 9)
#Figure 6 (page 10)
#Figure 7 (page 11)
#Figure 8 (page 12)
#Table A1 and A2 (Appendix)

      #Load Packages

install.packages("tidyverse")
library (tidyverse)

install.packages ("xtable")
library(xtable)

install.packages("forcats")
library(forcats)

install.packages ("data.table")
library(data.table)

install.packages ("ggthemes")
library(ggthemes)

install.packages("devtools")
library(devtools)

devtools::install_github("hrbrmstr/waffle")
library(waffle)

install.packages("RColorBrewer")
library(RColorBrewer)

install.packages ("cshapes")
library("cshapes")

install.packages ("sf")
library(sf)

install.packages("sp")
library(sp)

      #Import MMAD-RA Events Data, MMAD-RA Events Report Data and Original MMAD data

events <- read_csv("mmad_ra_events.csv")
eventreports <- read.csv("mmad_ra_eventreports.csv")
MMAD_events <- read_csv("mmadv1_events.csv")


                                  #Table 1 (page 5) - omitted

#create summary statistics for each actor type, if actor is 1
#police
pol <- events %>% group_by(police) %>% filter (police == 1) %>% summarise(mean(number_of_people_injured_num_mean, na.rm = T), sd(number_of_people_injured_num_mean, na.rm = T),
                                                 sum(is.na(number_of_people_injured_num_mean)), mean(number_of_people_killed_num_mean, na.rm = T),
                                                 sd(number_of_people_killed_num_mean, na.rm = T), sum(is.na(number_of_people_killed_num_mean)),  Count = n())
#milpolice
milpol <- events %>% group_by(milpolice) %>%  filter (milpolice == 1) %>%summarise(mean(number_of_people_injured_num_mean, na.rm = T), sd(number_of_people_injured_num_mean, na.rm = T),
                                                       sum(is.na(number_of_people_injured_num_mean)), mean(number_of_people_killed_num_mean, na.rm = T),
                                                       sd(number_of_people_killed_num_mean, na.rm = T), sum(is.na(number_of_people_killed_num_mean)),
                                                       Count = n())                                                
#military
mily <- events %>% group_by(military) %>% filter (military == 1) %>% summarise(mean(number_of_people_injured_num_mean, na.rm = T), sd(number_of_people_injured_num_mean, na.rm = T),
                                                    sum(is.na(number_of_people_injured_num_mean)), mean(number_of_people_killed_num_mean, na.rm = T),
                                                    sd(number_of_people_killed_num_mean, na.rm = T), sum(is.na(number_of_people_killed_num_mean)),
                                                    Count = n())
                                                    
#militia
milt <- events %>% group_by(militia) %>% filter (militia == 1) %>% summarise(mean(number_of_people_injured_num_mean, na.rm = T), sd(number_of_people_injured_num_mean, na.rm = T),
                                                   sum(is.na(number_of_people_injured_num_mean)), mean(number_of_people_killed_num_mean, na.rm = T),
                                                   sd(number_of_people_killed_num_mean, na.rm = T), sum(is.na(number_of_people_killed_num_mean)),
                                                   Count = n())

#sec forc
secf <- events %>% group_by(secforc_amb) %>% filter (secforc_amb == 1) %>% summarise(mean(number_of_people_injured_num_mean, na.rm = T), sd(number_of_people_injured_num_mean, na.rm = T),
                                                       sum(is.na(number_of_people_injured_num_mean)), mean(number_of_people_killed_num_mean, na.rm = T),
                                                       sd(number_of_people_killed_num_mean, na.rm = T), sum(is.na(number_of_people_killed_num_mean)),
                                                       Count = n())

#other
oth <- events %>% group_by(other) %>% filter (other == 1) %>% summarise(mean(number_of_people_injured_num_mean, na.rm = T), sd(number_of_people_injured_num_mean, na.rm = T),
                                                sum(is.na(number_of_people_injured_num_mean)), mean(number_of_people_killed_num_mean, na.rm = T),
                                                sd(number_of_people_killed_num_mean, na.rm = T), sum(is.na(number_of_people_killed_num_mean)),
                                                Count = n())
#missing
missing <- events %>% group_by(actor_missing) %>% filter (actor_missing == 1) %>% summarise(mean(number_of_people_injured_num_mean, na.rm = T), sd(number_of_people_injured_num_mean, na.rm = T),
                                                             sum(is.na(number_of_people_injured_num_mean)), mean(number_of_people_killed_num_mean, na.rm = T),
                                                             sd(number_of_people_killed_num_mean, na.rm = T), sum(is.na(number_of_people_killed_num_mean)),
                                                             Count = n())
#total                                                                                                                       
total <- events %>% summarise(mean(number_of_people_injured_num_mean, na.rm = T), sd(number_of_people_injured_num_mean, na.rm = T),
                              sum(is.na(number_of_people_injured_num_mean)), mean(number_of_people_killed_num_mean, na.rm = T),
                              sd(number_of_people_killed_num_mean, na.rm = T), sum(is.na(number_of_people_killed_num_mean)),
                              Count = n())
#Bind together
pol$ra_category<- "Police"
milpol$ra_category<- "Militarized Police"
mily$ra_category<- "Military"
milt$ra_category<- "Militia"
secf$ra_category<- "Ambiguous Security Force"
oth$ra_category<- "Other"
missing$ra_category<- "Missing"
total$ra_category <- "Total"

pol <- pol %>% select (-police)
milpol <- milpol %>% select (-milpolice)
mily <- mily %>% select (-military)
milt <- milt %>% select (-militia)
secf <- secf %>% select (-secforc_amb)
oth <- oth %>% select (-other)
missing <- missing %>% select (-actor_missing)

all <- rbind(pol, milpol) 
all2 <- rbind(all, mily) 
all3 <- rbind(all2, milt) 
all4 <- rbind(all3, secf) 
all5 <- rbind(all4, oth) 
all6 <- rbind(all5, missing) 
all7 <- rbind(all6, total)

df_table1 <- all7 %>% select (-Count) %>% select (ra_category, everything())
df_table2 <- df_table1 %>%  rename(Mean_Injured = `mean(number_of_people_injured_num_mean, na.rm = T)`,
                                    SD_Injured = `sd(number_of_people_injured_num_mean, na.rm = T)`,
                                    NA_Injured = `sum(is.na(number_of_people_injured_num_mean))`,
                                    Mean_Killed = `mean(number_of_people_killed_num_mean, na.rm = T)`,
                                    SD_Killed = `sd(number_of_people_killed_num_mean, na.rm = T)`,
                                    NA_Killed = `sum(is.na(number_of_people_killed_num_mean))`,
                                    Actor_Type = ra_category)
#Table to Latex code:
#Small design edits conducted in Latex (Overleaf) itself
xtable(df_table2)



                          #Figure 3 (page 6)
actor_type <- c('Police', 'Militarized Police', 'Military', 'Militia','Amb. Sec. Force',  'Other', 'Missing')
freqc <- all6$Count
df <- data.frame(actor_type, freqc)

fig3 <- df %>%
  mutate(actor_type = fct_relevel(actor_type, 
                                  "Police", "Militarized Police", "Military", 
                                  "Militia","Amb. Sec. Force",  "Other", "Missing")) %>%
  ggplot( aes(x=actor_type, y=freqc)) +
  geom_bar(stat="identity", width=0.4) +
  ylab("Number of MMAD-RA events") + xlab ("") + theme_bw()


                          #Figure 4 (page 7)
#filter events data per actor type
mily <- events %>% filter (military == 1)
mily$ra_category<- "military"

milt <- events %>% filter (militia == 1)
milt$ra_category<- "militia"

milpol <- events %>% filter (milpolice == 1)
milpol$ra_category<- "mil. police"

pol <- events %>% filter (police == 1)
pol$ra_category<- "police"

secfor <- events %>% filter (secforc_amb == 1)
secfor$ra_category<- "amb. security force"

oth <- events %>% filter (other == 1)
oth$ra_category<- "other"

missing <- events %>% filter (actor_missing == 1)
missing$ra_category<- "missing"

total <- rbind(mily, milt) 
total2 <- rbind(total, milpol) 
total3 <- rbind(total2, pol) 
total4 <- rbind(total3, secfor) 
total5 <- rbind(total4, oth) 
total6 <- rbind(total5, missing) 

# Omit NA's on injured and killed mean data

df_ra <- total6[!is.na(total6$number_of_people_injured_num_mean), ]                
df_ra2 <- total6[!is.na(total6$number_of_people_killed_num_mean), ]

#group data by actor type for injured and killed data
plotdata_inj <-  df_ra  %>%
  group_by(ra_category) %>%
  summarize(n = n(),
            median = median(number_of_people_injured_num_mean),
            Q1=quantile(number_of_people_injured_num_mean,probs = 0.25),
            Q3=quantile(number_of_people_injured_num_mean,probs = 0.75),
            sd = sd(number_of_people_injured_num_mean),
            se = sd / sqrt(n),
            ci = qt(0.975, df = n - 1) * sd / sqrt(n))
      
plotdata_kil <-  df_ra2 %>%
  group_by(ra_category) %>%
  summarize(n = n(),
            median = median(number_of_people_killed_num_mean),
            Q1=quantile(number_of_people_killed_num_mean,probs = 0.25),
            Q3=quantile(number_of_people_killed_num_mean,probs = 0.75),
            sd = sd(number_of_people_killed_num_mean),
            se = sd / sqrt(n),
            ci = qt(0.975, df = n - 1) * sd / sqrt(n))
#prepare for merge
plotdata_inj$type <- "Injuries"
plotdata_kil$type <- "Fatalities"

merge_kilinj <- rbind(plotdata_inj, plotdata_kil)
fig4data <- merge_kilinj %>%
  mutate(ra_category = fct_relevel(ra_category, "other",
                                   "police", "mil. police","missing", "amb. security force",
                                   "military","militia")) 


#plot figure 4 
fig4 <- ggplot(fig4data, 
               aes(x = ra_category, 
                   y = median, 
                   colour = type)) +
  geom_point(position=position_dodge(width=0.5), size = 3) + scale_colour_manual(values = c("Injuries" = "grey",
                                                                                            "Fatalities" ="black"),  guide = guide_legend(reverse = TRUE))  + 
  labs(x = "Repressive Actor", y = "Median Number of Injuries/Fatalities per Event") + 
  scale_x_discrete(labels = c('Other', 'Police', 'Militarized Police','Missing','Amb. Sec. Force', 'Military','Militia')) + 
  scale_y_continuous(breaks=seq(0,80,5)) + 
  geom_errorbar(aes(ymin = Q1, 
                    ymax = Q3), 
                width = .1, position=position_dodge(width=0.5)) + 
  theme_bw() + coord_flip() 
fig4 <- fig4 + theme(legend.title = element_blank())

#ggplot2::ggsave(filename= "CI_RA2_1102.pdf",
 #               plot = fig4,
 #               device = "pdf",
  #              width = 15,
   #             height = 15,
    #            units = "cm",
     #           dpi = 400)



                              #Figure 5a (page 9)
#prepare data
total6$date <- as.Date(total6$event_date)
fig5a_data <- total6 %>%
  mutate(ra_category = fct_relevel(ra_category, 
                                   "missing", "other","amb. security force", "militia", 
                                   "military", "mil. police", "police"))
#use original MMAD Events data here to plot the black line
MMAD_events$date <- as.Date(MMAD_events$event_date)

#plot Figure 5a
Fig5a <-  ggplot() + geom_freqpoly(data = MMAD_events, aes(x= lubridate::floor_date(date, "3 months")),binwidth = 90)  + 
  geom_bar(data = fig5a_data, aes (x = lubridate::floor_date(date, "3 months"), fill = as.factor(ra_category))) + scale_x_date(breaks = scales::breaks_pretty(10)) +  
  labs(x = "Time (quarter)", y = "Number of MMAD-RA events",
       fill="Repressive Actor") +  scale_fill_manual(values = c("police" = "#4E79A7", "mil. police" = "#F28E2B","military" = "#E15759","militia" = "#76B7B2",  "amb. security force" = "#59A14F","other" = "#EDC948","missing" = "#B07AA1"),  labels = c("Police", "Militarized Police", "Military", "Militia","Amb. Sec. Force", "Other", "Missing")) + 
  theme_bw()

                            #Figure 5b (Page 9)
#prepare data

stack <- total6 %>% group_by((lubridate::floor_date(date, "3 months")), ra_category) %>% summarise(count = n())
names(stack)[names(stack)== "(lubridate::floor_date(date, \"3 months\"))"] <- "quarterly"


fig5b_data <- stack %>% group_by(quarterly) %>% mutate(sum = sum(count)) %>% 
  group_by(ra_category, quarterly) %>% mutate(percentage = count/sum) %>% 
  select(-sum)


#plot Figure 5b
Fig5b <- ggplot(fig5b_data, aes(fill = factor(ra_category, levels=c("missing","other","amb. security force", "militia", "military","mil. police","police")), y=percentage, x=quarterly)) + 
  geom_bar(position="stack", stat="identity") +  scale_x_date(breaks = scales::breaks_pretty(10)) +
  labs(x = "Time (quarter)", y = "Percentage of MMAD-RA events") + ylim(0,1) +
  guides(fill = guide_legend(title = "Repressive Actor")) +
  scale_fill_manual(values = c("police" = "#4E79A7", "mil. police" = "#F28E2B","military" = "#E15759","militia" = "#76B7B2",  "amb. security force" = "#59A14F","other" = "#EDC948","missing" = "#B07AA1"),  labels = c("Police", "Militarized Police", "Military", "Militia","Amb. Sec. Force", "Other", "Missing")) + 
  theme_bw()

                          #Figure 6 (page 10)


# country-level data
country_eng <-  events %>%
  group_by(state_name) %>%
  summarize(n = n(),
            mean_engagement = mean(sec_engagement_MMAD))
          
#Use previous actor type data

RA_country <- data.frame(rbind(table(total6$state_name, total6$ra_category)))
RA_country2<- cbind(rownames(RA_country), data.frame(RA_country, row.names=NULL))
names(RA_country2)[names(RA_country2)=="rownames(RA_country)"] <- "state_name"

country <- left_join(RA_country2, country_eng, by= "state_name")

#Change n to include every ra observation

country$count <- (country$military + country$militia + country$mil..police + country$other + country$police + country$missing + country$amb..security.force)

country$pp_police_total <- country$police/country$count
country$pp_military_total <- country$military/country$count
country$pp_mp_total <- country$mil..police/country$count
country$pp_sf_total <- country$amb..security.force/country$count
country$pp_other_total <- country$other/country$count
country$pp_missing_total <- country$missing/country$count
country$pp_militia_total <- country$militia/country$count

#Change labeling of Côte d'Ivoire

country$state_name <- ifelse (country$state_name == "CÃ´te dâ€™Ivoire", "Côte d'Ivoire", country$state_name)

#long format
c21 <- country %>% filter (n >= 5) %>% select (state_name, mean_engagement, pp_police_total, pp_military_total, pp_other_total, pp_missing_total, pp_militia_total, pp_mp_total, pp_sf_total)
long2 <- melt(setDT(c21), id.vars = c("state_name"), variable.name = "ra_category")

#order the countries per level of engagement

agg2 <- aggregate(value ~ state_name, data=long2[long2$ra_category== "mean_engagement",], sum)
long2$state_name <- factor(long2$state_name, levels=agg2[order(agg2$value, decreasing = F), "state_name"])
long3 <- subset(long2, long2$ra_category != "mean_engagement")

#Reshuffle order
figure6_data <- long3 %>% 
  mutate(ra_category = fct_relevel(ra_category, 
                                   "pp_police_total", "pp_mp_total", "pp_military_total", 
                                   "pp_militia_total","pp_sf_total", "pp_other_total", "pp_missing_total"))%>%  as.data.frame()


#plot Figure 6

fig6 <- ggplot(figure6_data, aes(fill=ra_category, y=value, x=state_name, order = ra_category)) + 
  geom_bar(position="stack", stat="identity") + 
  labs(x = "Country", y = "Repressive actor distribution (percent)")  + ylim(0,1) +
  guides(fill = guide_legend(title = "Repressive Actor")) + coord_flip() + scale_fill_tableau(labels = c("Police", "Militarized Police", "Military", "Militia", "Amb. Sec. Force", "Other", "Missing")) +
  theme_minimal() 

                                                        #Figure 7 (page 11)
#Transform data to identify and get to the descriptions appearing in more than 20 events

events$id2 <- row.names(events)
ed3<- separate_rows(events, other_details_about_repressive_event, sep= ";")
ed3$event_details <- ed3$other_details_about_repressive_event
#tweak data to group the same description as the same description
ed3$event_details <- ifelse(ed3$event_details == "police used tear gas", "tear gas used", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == " tear gas used", "tear gas used", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == " arrests made", "arrests made", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == " stones thrown", "stones thrown", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == " property destroyed", "property destroyed", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == " shots fired", "shots fired", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == "rocks thrown", "stones thrown", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == "police used tear gas and water cannons", "tear gas used;water cannons", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == " water cannons", "water cannons", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == " tear gas", "tear gas used", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == " people arrested", "arrests made", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == "police made arrests", "arrests made", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == "police used batons", "batons used", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == "water cannons used", "water cannons", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == " police used tear gas", "tear gas used", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == "and used water cannons", "water cannons", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == "and water cannons", "water cannons", ed3$event_details)
ed3$event_details <- ifelse(ed3$event_details == " batons used", "batons used", ed3$event_details)

df3 <- ed3 %>% filter(event_details == "arrests made"|event_details == "tear gas used"|event_details == "shots fired"|event_details == "property destroyed"|event_details == "stones thrown"|event_details == "batons used"|event_details == "water cannons")

#now only filter unique id and event details combinations
#to not have two tear gas used for the same event etc.

events2 <- df3 %>% distinct(event_details, id2, .keep_all = "TRUE")
#now we have a dataset with the 7 identified event details per event 
#count per actor type to make waffles
#as such, first create actor type subsets such as done before

#subset and calculate frequency of string variable

mily2 <- events2 %>% filter (military == 1)
mily2$ra_category<- "military"
mily3 <- mily2 %>% select(event_details)
mily4 <- as.data.frame(table(mily3))

milt2 <- events2 %>% filter (militia == 1)
milt2$ra_category<- "militia"
milt3 <- milt2 %>% select(event_details)
milt4 <- as.data.frame(table(milt3))

milpol2 <- events2 %>% filter (milpolice == 1)
milpol2$ra_category<- "mil. police"
milpol3 <- milpol2 %>% select(event_details)
milpol4 <- as.data.frame(table(milpol3))

pol2 <- events2 %>% filter (police == 1)
pol2$ra_category<- "police"
pol3 <- pol2 %>% select(event_details)
pol4 <- as.data.frame(table(pol3))

secfor2 <- events2 %>% filter (secforc_amb == 1)
secfor2$ra_category<- "amb. security force"
sf3 <- secfor2 %>% select(event_details)
sf4 <- as.data.frame(table(sf3))

oth2 <- events2 %>% filter (other == 1)
oth2$ra_category<- "other"
oth3 <- oth2 %>% select(event_details)
oth4 <- as.data.frame(table(oth3))

mis2 <- events2 %>% filter (actor_missing == 1)
mis2$ra_category<- "actor missing"
mis3 <- mis2 %>% select(event_details)
mis4 <- as.data.frame(table(mis3))

#colors
myPalette <- brewer.pal(7, "Set2") 
names(myPalette) <- levels(pol4$event_details)
colScale <- scale_fill_manual(name = "",values = myPalette)


# transform data so its per 100 observations (every waffle 1 percent)

pol4$per <- pol4$Freq/sum(pol4$Freq)*100
pol4$Round_off <- round(pol4$per)

milpol4$per <- milpol4$Freq/sum(milpol4$Freq)*100
milpol4$Round_off <- round(milpol4$per)
#round shots fired down to 11, as the lowest value, to have 100 percent
milpol4$Round_off[4] <- 11

mily4$per <- mily4$Freq/sum(mily4$Freq)*100
mily4$Round_off <- round(mily4$per)

milt4$per <- milt4$Freq/sum(milt4$Freq)*100
milt4$Round_off <- round(milt4$per)

sf4$per <- sf4$Freq/sum(sf4$Freq)*100
sf4$Round_off <- round(sf4$per)
#round batons used down to 2, as the lowest value, to have 100 percent
sf4$Round_off[2] <- 2

mis4$per <- mis4$Freq/sum(mis4$Freq)*100
mis4$Round_off <- round(mis4$per)

oth4$per <- oth4$Freq/sum(oth4$Freq)*100
oth4$Round_off <- round(oth4$per)
#round tear gas up to 22 to have 100 percent
oth4$Round_off[6] <- 22

#Add everything together per category

m1 <- pol4 %>% select (event_details, Round_off)
m1$ra <- "Police"
m2 <- sf4 %>% select (event_details, Round_off)
m2$ra <- "Amb. Sec. Force"
m3 <- milpol4 %>% select (event_details, Round_off)
m3$ra <- "Mil. Police"
m4 <- mily4 %>% select (event_details, Round_off)
m4$ra <- "Military"
m5 <- milt4 %>% select (event_details, Round_off)
m5$ra <- "Militia"
m6 <- oth4 %>% select (event_details, Round_off)
m6$ra <- "Other"
m7 <- mis4 %>% select (event_details, Round_off)
m7$ra <- "Missing"

wf <- rbind(m1, m2) 
wf2 <- rbind(wf, m3) 
wf3 <- rbind(wf2, m4) 
wf4 <- rbind(wf3, m5) 
wf5 <- rbind(wf4, m6) 
wf6 <- rbind(wf5, m7) 

#reverse order
figure7_data <- wf6 %>%
  mutate(across(ra, factor, levels=c("Police",  "Mil. Police", "Military","Militia", "Amb. Sec. Force",
                                     "Other", "Missing")))
#plot Figure 7


waffles_fig7 <- ggplot(figure7_data, aes(fill=event_details, values=Round_off)) +
  geom_waffle(color = "white", size=1, n_rows = 10) +
  facet_wrap(~ra, ncol=3) + 
  scale_x_discrete(expand=c(0,0)) +
  scale_y_discrete(expand=c(0,0)) +
  scale_fill_tableau(name=NULL) +
  coord_equal() +
  theme_enhance_waffle(
  ) + theme_minimal() + theme(strip.text.x = element_text(
    size = 10,  face = "italic")) +  theme(legend.position="bottom")

                                                                                            ##Figure 8 (page 12)


#Use previous actor type - event data 
#filter country and year

MENA <- total6 %>% filter (cowcode == 600| cowcode == 615 | cowcode == 616 | cowcode == 620 | cowcode == 651 | cowcode == 652)
MENA$id <- rownames(MENA)
MENA$event_date <- as.Date(MENA$event_date, "%Y-%m-%d")
MENA <- MENA %>% filter (event_date > '2009-12-31')

#Rely on original MMAD data to connect location code to city name
names <- MMAD_events %>% filter (cowcode == 600| cowcode == 615 | cowcode == 616 | cowcode == 620 | cowcode == 651 | cowcode == 652) %>%  select(location, latitude, longitude, asciiname)

#merge to get the latitude/longitude and city names
names <- distinct(names)
names2 <- names %>% distinct (location, .keep_all = TRUE)

NAM3 <- left_join(MENA, names2, by = "location")

worldmap2010 <- cshp(date=as.Date("2010-06-30"), useGW=F)
border.mor <- worldmap2010[worldmap2010$cowcode==600 | worldmap2010$cowcode==615 | worldmap2010$cowcode==616 | worldmap2010$cowcode==620|worldmap2010$cowcode==651|worldmap2010$cowcode==652,]

#functions
calculate_mode <- function(x) {
  uniqx <- unique(na.omit(x))
  uniqx[which.max(tabulate(match(x, uniqx)))]
}

mean.new <- function(v) {
  if (all(is.na(v))) { return(NA) } else { return(mean(v, na.rm=T)) }
}

#location data
MD <- NAM3 %>% group_by(location, longitude, latitude) %>% summarize(number = n(), actor = calculate_mode(ra_category),
                                                                     mean_numinjured_perlocation = mean.new(number_of_people_injured_num_mean), mean_numkilled_perlocation = mean.new(number_of_people_killed_num_mean))
fig8_data <- MD %>% filter (!is.na(longitude)) %>%
  st_as_sf(coords = c("longitude", "latitude"), crs=4326)

#Plot Figure 8

cities <- names2 %>% filter(asciiname == "Cairo"| asciiname == "Tunis"| asciiname == "Benghazi"| asciiname == "Rabat"| asciiname == "Algiers" | asciiname == "Tripoli" | asciiname == "Damascus" | asciiname == "Aleppo")
y <- ggplot(border.mor) + geom_sf(fill = "white") 
y2 <- y + geom_sf(data = fig8_data, aes(size = number, color = actor)) + geom_text(data = cities, aes(x = longitude, y = latitude, label = asciiname), size=3, color = "black")
y3 <- y2 + labs (fill = "actor")
y4 <- y3 + theme(panel.grid.major = element_blank(), 
                 plot.margin=unit(c(0.5,0.5,0.5,0.5),"cm")) + guides(color = guide_legend(title = "Repressive Actor    ")) + 
  scale_color_manual(values = c( "police" = "#4E79A7", "mil. police" = "#F28E2B","military" = "#E15759","militia" = "#76B7B2",  "amb. security force" = "#59A14F","other" = "#EDC948","missing" = "#B07AA1"),  
                     labels = c("Police", "Mil. Police", "Military","Militia","Amb. Sec. Force", "Other", "Missing")) + guides(size = "none" )
fig8 <- y4 +
  theme(axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(),  
        axis.ticks.y=element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(), panel.background = element_rect(fill='white', colour='white'), legend.position="bottom")

                                                                                        ##Table A1 and A2 (Appendix)

#Load packages 
#put here to not mask other packages

install.packages("expss")
library(expss)

install.packages("vtable")
library(vtable)


#Event dataset
newdatalabels <-apply_labels(events,
                             event_date="Event date",
                             location ="City ID in the Geonames Database",
                             cowcode ="Correlates of War country numeric",
                             state_name = "Country name",
                             issue_MMAD = "Issue/motivation for protest",
                             actors_MMAD = "The actors involved in the protest",
                             side_MMAD= "The side towards which the protest is directed at",
                             scope_MMAD = "The government level towards which the protest is directed at",
                             part_violence_MMAD = "Level of violence from protest participants",
                             sec_engagement_MMAD = "Level of official security forces engagement",
                             repressive_actor_name = "The name of the repressive actor",
                             other_details_on_repressive_actor = "Other details on the repressive actor",
                             other_details_about_repressive_event = "Other details about repressive event",
                             milpolice  = "Repressive actor type is military police",
                             military  = "Repressive actor type is military",
                             militia  = "Repressive actor type is militia",
                             police = "Repressive actor type is police",
                             secforc_amb  = "Repressive actor type is ambiguous security force",
                             other  = "Repressive actor type is other",
                             number_of_people_injured  = "Number of people injured",
                             number_of_people_injured_num_min = "Number of people injured, low",
                             number_of_people_injured_num_mean = "Number of people injured, mean",
                             number_of_people_injured_num_max = "Number of people injured, high",
                             number_of_people_killed = "Number of people killed",
                             number_of_people_killed_num_min= "Number of people killed, low",
                             number_of_people_killed_num_mean  = "Number of people killed, mean",
                             number_of_people_killed_num_max   = "Number of people killed, high",
                             num_participants_MMAD  = "Number of participants",
                             num_participants_num_min = "Number of participants, low",
                             num_participants_num_mean = "Number of participants, mean",
                             num_participants_num_max = "Number of participants, high",
                             actor_missing = "Information on repressive actor type is missing",
                             nr_actors = "Number of actors")

vt(newdatalabels,out='latex',class=F, file='vartable.tex')
#Manually added column indicating source of the variable in Latex
#Repeat for event reports dataset

newdatalabels2 <-apply_labels(eventreports,
                              id = "Event report ID",
                              event_date="Event date",
                              location ="City ID in the Geonames Database",
                              cowcode ="Correlates of War country numeric",
                              state_name = "Country name",
                              issue_MMAD = "Issue/motivation for protest",
                              actors_MMAD = "The actors involved in the protest",
                              side_MMAD = "The side towards which the protest is directed at",
                              scope_MMAD = "The government level towards which the protest is directed at",
                              part_violence_MMAD = "Level of violence from protest participants",
                              sec_engagement_MMAD = "Level of official security forces engagement",
                              repressive_actor_name = "The name of the repressive actor(s)",
                              other_details_on_repressive_actor = "Other details on the repressive actor(s)",
                              other_details_about_repressive_event = "Other details about repressive event",
                              milpolice  = "Repressive actor type is military police",
                              military  = "Repressive actor type is military",
                              militia  = "Repressive actor type is militia",
                              police = "Repressive actor type is police",
                              secforc_amb  = "Repressive actor type is ambiguous security force ",
                              other  = "Repressive actor type is other",
                              number_of_people_injured  = "Number of people injured",
                              number_of_people_killed = "Number of people killed",
                              num_participants_MMAD  = "Number of participants",
                              actor_missing = "Information on repressive actor type is missing",
                              summary_of_event = "Summary of event") 

vt(newdatalabels2,out='latex',class=F, file='vartable2.tex')
#Manually added column indicating source of the variable in Latex